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Abstract 

An implicit method for the computation of unsteady flows on unstructured grids 
is presented. Following a finite difference approximation for the time derivative, the 
resulting nonlinear system of equations is solved at each time step by using an agglom- 
eration multigrid procedure. The method allows for arbitrarily large time steps and is 
efficient in terms of computational effort and storage. Inviscid and viscous unsteady 
ffows are computed to validate the procedure. The issue of the mass matrix which 
arises with vertex-centered finite volume schemes is addressed. The present formula- 
tion allows the mass matrix to be inverted indirectly. A mesh point movement and 
reconnection procedure is described that allows the grids to evolve with the motion of 
bodies. As an example of flow over bodies in relative motion, flow over a multi-element 
airfoil system undergoing deployment is computed. 
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1 Introduction 


Solution techniques for computing steady flows on unstructured grids have evolved to a high 
degree of sophistication. With explicit schemes, convergence to steady state is usually un- 
acceptably slow, especially as the problem sizes and complexities grow. Therefore, either 
multigrid methods [20, 29] or implicit schemes [36, 44, 2, 4] are required to accelerate the 
convergence. On the other hand, solution techniques for dealing with unsteady flows have 
lagged behind. Explicit schemes, deemed to be too slow for obtaining steady state solutions, 
may be the schemes of choice for certain unsteady applications, when the time scales of 
interest are small, or more precisely, when they are comparable to the spatial scales. The 
grids should be clustered only in regions of interest; otherwise, the size of the explicit time 
step could become unnecessarily small. However, when dealing with many low frequency 
phenomena such as flutter, explicit schemes lead to large computing times. A time-accurate 
residual averaging formulation [18, 42], and temporal adaptation techniques [13], which en- 
able different cells to take a varying number of local time steps to get to a particular time 
level, can be used to realize modest improvements in the performance of explicit methods, 
but the sizes of the time steps are still controlled by the spatial resolution. Therefore, it is 
desirable to develop a fully implicit method, where the time step is solely determined by the 
flow physics and is not limited by the cell sizes. Also, for many practical viscous flows, the 
time step restrictions imposed by small cells deep inside the boundary layer are excessively 
small. Since the boundary layer is quasi-steady, implicit methods that allow for larger time 
steps may be more suitable. 

When an implicit scheme is used to compute unsteady flows, one has to drive the unsteady 
residual to zero (or at least to truncation error) at each time step. In the context of factored 
implicit schemes, this is usually done by employing inner iterations [31, 30, 7, 33]. It is the 
role of these inner iterations to eliminate errors, if any, due to factorization and linearization, 
and sometimes also errors arising from employing a lower order approximation on the implicit 
side. The number of inner iterations required may be large depending on the flow situation, 
the size of the time step employed and the extent of mismatch of the explicit and implicit 
operators. Jameson [11] has advocated the use of a full approximation storage multigrid 
procedure as a driver for a fully implicit scheme when using structured grids. The significant 
advantage of the approach when multigrid is used to solve the nonlinear problem is that 
it incurs no storage overheads plaguing traditional implicit schemes based on linearization. 
The method is therefore particularly attractive for unstructured grid computations in three 
dimensions. The method allows the time step to be determined solely based on flow physics. 
It has been used to compute two- and three-dimensional inviscid flows over airfoils and wings 
[11, 1] using structured grids. Vassberg [41] has applied this technique to compute flow 
solutions over oscillating airfoils using unstructured grids where a sequence of triangulations 
was generated by removing points from the fine grid triangulation. 

Multigrid techniques have been successfully extended to deal with unstructured grids 
by generating a sequence of non-nested grids and using piecewise-linear transfer operators 
[20, 29]. An important development in the use of multigrid techniques for unstructured grids 
is the agglomeration multigrid algorithm [15, 37, 43] for the two- and three-dimensional 
Euler equations. This technique has since been extended to deal with viscous flows [14, 22], 
The main advantage of the agglomeration multigrid algorithm is that it only requires the 
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fine mesh to be generated; the coarse grids are automatically generated by fusing fine grid 
control volumes using an efficient graph-based algorithm, resulting in a fully nested sequence 
of coarse grid levels. 

In this work, the agglomeration multigrid strategy is used to solve the nonlinear system 
of equations at each time step. The nested property of the agglomeration approach enables 
a straight-forward treatment of problems involving moving meshes. It is also shown that the 
mass matrix can be inverted indirectly during the multigrid process. The issue of the mass 
matrix is critically examined in both one and two dimensions. In order to allow the grids to 
conform to moving geometries, a technique is proposed and tested that attempts to maintain 
the validity and quality of the triangulation. Inviscid flow over a pitching airfoil and viscous 
flow over an impulsively started cylinder are computed and the results are compared with 
experiments and other computations. Finally, an exploratory computation is carried out 
that simulates the phenomenon of flap deployment. 


2 Governing equations and discretization 

The equations governing compressible fluid flow in integral form for a control volume V(t) 
with boundary S(t) are given by 

/ Wdv + i [F(W, n, s) — G(W, V1F, n)]da = 0, (1) 

Ot Jv(t) Js(t) 

W = [p,pV,pe] T 
F(W, n, s) = (V - s).nW, 

G(W, V1F, n) = [0, t,t-V - q.nf, 

where t and q are respectively, the stress and heat flux vectors, p is the density, V is the 
velocity vector with Cartesian components Vy e is the specific total energy, n is the outward 
unit normal vector of the boundary S(t) } and s is the velocity vector of the boundary. The 
equations are augmented by the equation of state, which for a perfect gas is 

p\V \ 2 

P=(l-l)(pe- P -^). (2) 

In the present scheme, the variables are stored at the vertices of a mesh composed of 
triangles. The control volumes are nonoverlapping polygons which surround the vertices of 
the mesh. The contour integrals in Eqn. (1) are replaced by discrete path integrals over the 
faces of the control volume which are computed using the trapezoidal rule. This technique 
can be shown to be equivalent to using a piecewise-linear finite-element discretization under 
certain conditions. For dissipative terms, a blend of Laplacian and biharmonic operators 
is employed [20]. The Laplacian term acts only in the vicinity of shocks and is inactive 
elsewhere, while the biharmonic term acts only in regions of smooth flow. Only the Laplacian 
dissipative term is used on the coarse grids. 

After applying the finite volume procedure, the following system of coupled differential 
equations is obtained: 

^-(VMW) + R(W) = 0 (3) 
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Here W is the solution vector over the whole held, RiW) is the residual vector approximating 
the second integral in Eqn. (1), V is the cell volume associated with the vertex and M is the 
mass matrix. 

The mass matrix arises because the update indicated by the residual RiW) should be 
made to the average value in the control volume. It thus relates the average value of a 
control volume associated with a vertex to the point values of the vertex and those of its 
immediate neighbors. This definition differs from the way the mass matrix is defined in finite 
element formulations, where the mass matrix arises naturally from requiring the residual of 
the discretized PDE to be orthogonal to a set of trial functions, with the solution expanded in 
a set of basis functions. If the dissipative terms are made proportional to the residuals, this 
definition of the mass matrix carries through, e.g. the Streamline Upwind Petrov Galerkin 
method [10]. For finite volume schemes employing a polynomial reconstruction procedure 
within a cell, we instead derive the mass matrix entries by computing the average of this 
polynomial over the control volume. The mass matrix M couples the system of ODE’s in 
Eqn. (3). The effect is that even with an explicit scheme, one has to deal with the solution 
of a coupled linear system. A technique called “mass-lumping” [38], replaces the matrix M 
by the identity matrix. For second-order accurate cell-centered schemes, which employ the 
triangles as the control volumes and store the values at the centroids, mass-lumping does not 
compromise the accuracy, since the point value at the centroid matches the average value 
to second order. However, for cell-vertex schemes on nonuniform grids, the centroid of the 
control volume is not represented by the vertex in question. For time-accurate computations 
on such grids, mass lumping would appear to introduce locally a first order spatial error. 
This approximation is routinely adopted for unsteady flows as well, and does not appear to 
adversely affect the quality of the solutions obtained. Davis and Bendiksen [9] observed little 
discernible differences in the unsteady solutions when using the full and the lumped mass 
matrices. However, since they used an explicit scheme, the time steps were quite small and 
furthermore, the grids employed appeared to be fairly uniform. The technique employed to 
solve the mass matrix (a few Jacobi iterations) in [19, 9] is not efficient, especially when larger 
grids are used. However, Miller [24] and Wathen [45] have established mesh-independent 
bounds on the eigenvalues of the diagonally preconditioned mass matrix arising out of linear 
finite elements and have reported that a conjugate gradient method can be used to efficiently 
invert the diagonally preconditioned mass matrix in just a few iterations. When higher order 
spatial discretizations are employed, the mass matrix has to be reckoned with, even when 
using cell-centered discretizations. We note that the mass matrix can be avoided altogether 
if only cell averages are employed for the spatial discretization. 


3 The implicit scheme 

We first outline the implicit scheme as developed by Jameson [11] for cell-centered, structured 
grids, where mass lumping was used. Replacing the mass matrix in Eqn. (3) by the identity 
matrix and making a 3-point backward-difference approximation for the time derivative 
yields 


V n+1 W n+1 

2 At 


V n W n + 
At 


y n - 1 IW “ 1 

2 At 
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+R(W n+1 ) = 0. 


( 4 ) 


When applied to a linear differential equation of the form, 


dW 

dt 


aW 


( 5 ) 


this discretization is A-stable i.e., stable for all values of a At in the left-half of the complex 
plane [8]. Eqn. (4) is now treated as a steady state equation by introducing a pseudo-time 
variable t*. The multigrid scheme then solves the following nonlinear system to steady state 
using local time steps At*: 

dVU 

— + nr( u ) = o, ( 6 ) 

where U is the approximation to W n+1 . Here the unsteady residual R*{U ) is defined as 

R*(U) = ±VU + R(U) - S{V n W n , W 1 ^- 1 ) (7) 


with the source term 


S{v n w n , y^r -1 ) 


— V n w n — f/ n - 1 H / n_1 

At 2 At 


( 8 ) 


remaining fixed through the multigrid procedure. 

A multi-stage Runge-Kutta scheme applied to solve Eqn. (6) performs the role of a 
smoother. A low-storage second order accurate m-stage Runge-Kutta scheme to advance U 
is given by 


Qo = U l 

V n+1 Q k = V n+1 Q 0 -a k At*R*(Q k _ 1 ) (9) 

U l+1 = Q m 

Starting with U 1 = W n , the sequence of iterates U l ,l = 1,2,3.... converges to W n+1 . 

However, the way the scheme has been formulated has been observed by Arnone et al. [3] 
to be unstable for small physical time steps, At. This is counter-intuitive because when 
using a small At, the multigrid procedure should converge fast and ideally, in the limit of 
explicit time steps, the multigrid procedure should converge in just a few iterations. Melson 
et al. [23] showed that the problem is due an instability that arises when a small At is used. 
They modified the scheme to get rid of this instability. The source of the problem is that 
the unsteady residual R*(W) includes the term and, is therefore, treated explicitly 

in the Runge-Kutta scheme. Their analysis showed that if this term were treated implicitly 
in the Runge-Kutta scheme, the stability region would grow as At is decreased. It is easy to 
treat the term implicitly since it is only a diagonal term. Splitting the residual R*{U) as 

r-(U) = 1 ^ t vu + R(U)-s, (io) 
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the Runge-Kutta scheme now becomes, 

Qo = U l 

I + ^a*Ar] v n+1 Qk = V n+1 Qo (11) 

- a k At* [R(Qk-i) — S'] 

U l+1 = Q m 

With the modihed scheme, Melson et al. [23] have shown that arbitrarily large or small At 
may be employed. 

As in [11, 23], we employ a full approximation storage multigrid scheme. The source 
term is computed only on the fine grid and the coarse grid problems are driven by the fine 
grid residuals. For the generation of coarse grids we follow the agglomeration multigrid 
procedure. In this method, a sequence of coarse grids is generated using efficient graph- 
based algorithms. This method has certain advantages when dealing with rigidly moving or 
deforming meshes. Since the edges that comprise the coarse grid volumes are subsets of the 
fine grid control volume edges, when the grid moves rigidly or deforms, the projections of 
the control volume faces onto the coordinate directions are easily computed from those of 
the fine grid. Also, as long as no grid points are added or removed, and the triangulation 
remains valid, and the grid connectivity remains unchanged, the interpolation operators 
are unchanged. Multigrid schemes based on non-nested triangulations would require the 
recomputation of the interpolation operators when the grids deform. Even if the mesh 
connectivity changes, the agglomeration algorithm is efficient enough that it can be used to 
regenerate the coarse grids without incurring substantial overheads. 

4 Treatment of the mass matrix 

When employing a vertex-centered approximation, making a 3-point backward-difference 
approximation for the time derivative yields 

3 2 

V n+1 M n+1 W n+1 V n M n W n 

2 At At 

+ -^ t V n - 1 M n ~ 1 W n - 1 + R(W n+1 ) = 0. (12) 

The multigrid scheme now solves the following system to steady state using local time steps 
At*: 

dVlJ 

— + nr(u) = o, (is) 

where U is the approximation to W n+1 , and R*(W) now includes the mass matrix terms. 
Notice that the first term ~§^VU does not involve the mass matrix, thus uncoupling the 
system of equations. The explicit Runge-Kutta scheme can be applied exactly as before. 
The inversion of the mass matrix is thus accomplished indirectly during the multigrid pro- 
cedure. However, the modihed scheme of Melson et al. [23] poses a serious problem. Their 
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modification would require the term j^VMU, which is no longer a diagonal term, to be 
treated implicitly. We have devised a modification that solves this problem which is detailed 
below. The implicit Runge-Kutta scheme that is stable for all At is given by 

Qo = U l 


1 + 


3 

2A t 


a k At*M n+1 


V n+1 Q k = V n+1 Q 0 - 

a k At* [R(Q k _i) — S'] 


U l+1 = Q 


(14) 


where the source term S is now given by 


S 


2 1 

— V n M n W n y n - 1 AP'“ 1 lT™'“ 1 

At 2At 


( 15 ) 


If we simply replace the mass matrix M by the identity on the left hand side of Eqn. (14), we 
have observed that the instability at small time steps persists. In our modification, we first 
add and subtract ^^a k At* M n+1 V n+1 Q k -\ on the right hand side of Eqn. (14) to obtain 


1 + 


2A t 


a k At*M 


* M n + 1 


V n+1 Q k = V n+1 Q 0 - 


akAt'R'iQk-!) + —a k At*M n+1 V n+1 Q k . 1} 


where use has been made of the equation 


R*{U) = 7^ t VMU + R ( U ) ~ S ■ 


( 16 ) 


(U) 


Note that the same term ^^a k At*MVQ appears on the left and the right hand sides of 
Eqn. (16), except that they are evaluated at the k — 1 and k stages, and that R*{U ) is being 
driven to zero. The mass matrix M can now be replaced by /3/, where / is the identity 
matrix and f3 is a constant yielding the following equation: 


1 + 


2A t 


a k At*/3 


V n+1 Q k = V n+1 Q C 


-a k At*R*(Q k - 1 ) + — akAt'pV'+'Qk-! 


( 18 ) 


The method can always be stabilized by increasing f3 and is akin to using a damped Jacobi 
method. The Runge-Kutta scheme no longer requires a matrix inversion. For small time 
steps of the order permitted by the explicit scheme, we find that the choice of f3 = 2 stabilizes 
the scheme. 
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5 Mesh point movement strategies 

In order to be able to perform unsteady flow simulations over moving geometries, a body- 
conforming mesh has to be regenerated either globally at each time step, or the existing grid 
can be allowed to deform. The former option is expensive, especially in three dimensions. 
However, Baum et al. [6] have proposed some simplifying strategies. Whenever the body 
movement becomes too severe, they regenerate a coarse mesh either locally or globally. This 
is followed by the use of adaptive h-rehnement in order to create a fine mesh. In the present 
work, we only investigate mesh point movement strategies. Tension spring analogy [5, 28, 35] 
or other physical analogies, such as incompressible flow [12], are typically used to move the 
mesh points. In the former case, the distribution of the spring stiffnesses is crucial. Since the 
techniques try to maintain the connectivity of the grid at all costs, the grid lines may cross 
resulting in invalid triangulations. Nevertheless, for many simple configurations, such as 
isolated airfoils, no crossover occurs and the spring analogy has been demonstrated to work 
well e.g., [5]. For such cases however, the use of exponentially varying scaling factors [9] or 
rigid body motion are simpler. The real challenge for any mesh point movement strategy is 
when relative motion is present between bodies in close proximity and when fine grids are 
employed. Other possibilities for grid movement are methods in use in the moving finite 
element method [25], where evolution equations are derived for grid point motion from the 
governing equations. 

In our approach, we use the tension spring analogy to allow the grid points to react to 
the motion of the geometries. The following linear system of equations is solved by a Jacobi 
method: 

X j) — Siji (19) 

j 

where the source term Sij is computed so as to maintain the initial grid in the absence of any 
displacements. The spring stiffness Kij is taken to be !- p where is the length of the edge 
joining nodes i and j. Over a number of applications, we have found the value of p = 2 to 
work well. When the boundary motion is prescribed, this method does not guarantee that 
the grid lines will not cross. The next improvement is to make the spring system nonlinear 
i.e. , the boundary motion is decomposed into smaller steps and the procedure is repeated at 
every step. When relative motion is present, this results in excessive skewing of grid lines and 
eventually the lines do cross. Thus some kind of reconnection procedure is unavoidable. We 
choose to swap the edges by the Delaunay criterion. However, even though the grid remains 
valid, the resulting point distribution is typically unsatisfactory. Therefore a presmoothing 
procedure is used to smooth the distribution of points. This involves using a Jacobi method 
on the following system: 


(. I + eNi)x? ew = xf + e]T x 3 new } (20) 

j 

where W is the degree of node i. Typically we find that of the order of 100-200 steps 
are required to solve Eqn. (19), while only about 3-4 iterations of Eqn. (20) are sufficient. 
The number of Jacobi iterations to solve Eqn. (19) may appear to be large, but this is 
mainly due to the large displacements arising from using large time steps permitted by the 
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implicit scheme. The grid motion procedure described above also has applications in design 
optimization, where the surface geometry changes during the design cycle. 

For the implicit flow solver, the coarse grids are rederived by using the agglomeration 
algorithm whenever the grid connectivity changes. An incremental agglomeration algorithm 
is possible that operates only on the affected region, but this is not currently done. Typically, 
for the time steps used with the implicit scheme, the grid reconnection and reagglomeration 
consume about a third of the time required for the flow solver, mainly on account of the 
swapping procedure not being vectorized. 

The grid movement terms need to be discretized carefully so that freestream is preserved. 
In other words, simply moving the grid through the domain should not change the freestream 
solution. The Geometric Conservation Law (GCL) [39, 47, 27] formalizes this concept. It 
can be derived from the continuity equation in Eqn. (1) by first assuming the control volumes 
to be the simplices themselves. Assuming a uniform velocity held and a constant density 
held, we obtain 

^ + i M [v - 8] - nda = 0 ' (21) 

where V is the velocity held and s is the velocity of the boundary S(t). Since V is constant 
and the control volume is assumed to be closed at all times so that n da = 0, the 
equation becomes 

dV f 

— f s.n da = 0. (22) 

dt Js(t) v ; 

The discrete form of this equation should hold at all time steps and for all the simplices and 
is called the GCL. Using a forward Euler approximation for the time derivative, we obtain 


ya+l-ya = 


rt+At 


s.n da dt 


= e r ' 


s.n da dt . 




( 23 ) 


where Si(t) = Yj U/j(t) i s the surface enclosing the volume Vi(t) of simplex /. The term 
inside the summation represents the volume swept out by the boundary fi/j as the grid 
points forming that segment move. If the grid points are allowed to move arbitrarily, the 
GCL enables the edge velocity s and the grid normal n to be determined. Since simplices 
are convex, the volumes V n , V n+1 are uniquely determined by the positions of the points at 
time levels n and n + 1. In two dimensions, it can be shown [47] that the change in volume 
can be expressed as 

V n+i _ V n = At j2 s n 3 +1/2 .N] +1/2 , (24) 

3 

where the summation is over the edges forming the triangle I. Here the edge velocity s™ +1 ^ 2 
is given by the average of the velocities of the vertices connected by the edge j, and IVj, the 
normal scaled by the edge length, is given by the average of the values at the old and new 
time level: 

N] +1/2 = 0.5 [N] + N] +1 ). (25) 
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The velocity is assumed to be Constant within a time step. Therefore, the velocity of vertex 
i is given by 


n+l/2 


y ; >+1 - 

Ai 



(26) 


where X is the position vector of the vertex. 


j 



Figure 1: Control volume for vertex i and edge normals. 


For a vertex-centered scheme, such as the one used in the present work, the control 
volumes are formed by l lie median dual. We closely follow the excellent development of 
Nkonga and Guillard [27], who have presented the steps involved in properly discretizing the 
flow equations in three dimensions so as to obey the GCL for a two-level scheme in time. The 
control volume edges are delimited by segments of the medians of the triangles (Figure 1). 
Inspecting one fedge of the triangle, we notice that the control volume edge is comprised of 
two, generally non-collinear median segments. The results from Eqns.(24, 25) can be used to 
express the change in the volume of any polygon as a summation over the edges comprising 
the polygon of the volumes swept out. We obtain the following expression for the change in 
Vi, the volume associated with vertex i: 


T/ra+l T rn Ai W 7 yn+l/2 ra+1/2 . j^j'n+1/2 ra+1/2 

G — G - A ' 2^ •'../• •* ./. “T lyl j,R - S j,R ■> 

j 


(27) 


where the summation is over the edges meeting at vertex i. Also, N and the 

normals scaled by the lengths to the left and right of the edge, are obtained as averages of 
the values at n and n . + 1 time levels. The edge velocities, s J^ 2 and are computed as 

the averages of the velocity of the midpoint of the edge and that of the left or right centroid, 
respectively. 

In a vertex-centered finite- volume setting, the discretization of the governing equations 
(Ecpi. (1)) for a two-level explicit or implicit scheme that obeys the GCL is performed as 
follows: 


(MVW ) n+1 - ( MVW) n 
At 


+ 

+ 


]riw [n ]+ 1/2 .{v v - s n 3 i 1/2 ) 

j 

»,(»'". N" +1/2 ) = 0, 


+ i\c+ i/2 .(v - y,i 1/2 ( 


(28) 
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where rj = n for forward and r] = n + f for backward Euler methods, and W is the average of 
the values at the vertices connected by the edge j . Also, g 3 represents the discretization of 
the last term in Eqn. (1) and makes use of the scaled normal for the edge, N 3 +1 ^ 2 , defined 
as 

N] +1/2 = N n 3 l 1/2 + iVj+ 1/2 . (29) 

When a three point backward difference approximation is used for the time derivative, 
the spatial discretization is performed as follows: 

3(MVW) n+1 - A(MVW) n + {MVWy- 1 

2At 

+ W n+1 {-[N]+ L 1/2 .(V n+1 - 4t 1/2 ) + N]+ R 1/2 .(V n+1 - 4+ 1/2 )] 

3 

- )[Jv;h /2 w" +1 - yi 1/2 ) + n ]~ r i , 2 .( v '+' - s ;y /2 )]} 

+ ^ ) (»r” + rAT" +1/2 )-t 9) (»™+ 1 .Ar"- 1/2 ) = o. (30) 

It is easy to see that this formulation satisfies the GCL, and that it reduces to the standard 
three-point formula (Eqn. (12)) in the absence of any grid motion. We mention here that 
for fixed grid applications it is not necessary to compute N j,L and ^ j,R separately. Instead 
only the sum of the two is required, which can be computed easily from the the centroidal 
coordinates. On the other hand, for moving grid problems, the N j,L and ^ j,R need to be 
computed separately for use in Eqns. (28,30). 

The spatial discretization given above only deals with the convective and diffusive fluxes. 
The dissipative fluxes are computed in the usual fashion using the scaled normal given by 
Eqn. (29) for the two-level explicit or implicit scheme. For the three-point scheme, the scaled 
normal is defined as follows: 

Nj = 3/2 [lVj+ 1/2 + 1VJ+ 1/2 ] - 1/2 [N]~ L 1/2 + N]~ r 1/2 ] . (31) 

The discretization of the dissipative fluxes has no implications for the GCL as these fluxes 
vanish for a uniform held. 

We have verified that the formulation given above preserves the freestream conditions 
to machine precision. With other treatments, such as simply modifying the fluxes at the 
vertices by using nodal velocities, freestream conditions are preserved to much less precision, 
which could be detrimental in some applications. 

Finally, we address the issue of boundary conditions for inviscid hows. The control volume 
for a boundary vertex is illustrated in Figure 2. The summation over the edges in Eqns. 
(28,30) is augmented with the contributions from the the boundary edges that delimit the 
control volume. For the two-level scheme, the following terms are added to Eqn. (28): 

wi\Niy r ‘.(vi - yy /2 )] + w'aj\yy /2 .(n - <g /2 )] 

+ 9(W' L ",JVy 1/2 ) + <,(H-2,JVy 1/2 ), (32) 

where Wl and Wr are the values of W on the left and right side of vertex z, IV/,. r and IVyL 
are the outward normals to the boundary edges of the control volumes, and and Sy# 
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are the grid velocities associated with these edges. For an inviscid all, since the [V — s] = 0, 
only the last two terms are retained. Likewise for the three-point scheme (Ecpi. (30)), the 
additional terms are 


W£ +1 


0<i 1/2 .( v u l *' - <1 1/2 ) - - yyy 


+ IKr ,2 '( v r' - <U 2 ) - 5Jv;,; 1/2 .n 
+ fw » 2 +1 , n ;*' 12 ) + g ( w ?' , i w ;* L ' 12 )] 


‘•(vr-<k /2 ) 


1 + 1 


- y 9 (H2 +1 , jv"y /2 ) + 9(1-1^+', jv;y' 2 )] 


n+ 1 l\T n ~^-/' 2 ) 


For an inviscid wall, only the last four terms are retained, implying the following boundary 
condition: 


T-n+l/2 7i-|-l/2 


jn— 1/2 n— 1/2 


6 Mesh restructuring and interpolation 

The procedure for mesh restructuring was discussed earlier. This involved a mechanism for 
moving grid points in response to the motion of the boundary points, evolving the solution 
with the grid movement terms present, and swapping the edges of the triangulation to 
improve grid quality. When using the three point difference formula (Eqn.( 12) ), the swapping 
of edges at a particular time step has to be done in such a manner that the triangulation is also 
valid (i.e., no crossing of edges) at the previous time step. This additional criterion can easily 
be incorporated in the swapping procedure. Using the trapezoidal rule would circumvent 
this step and still ensure second order accuracy in time. However it is unattractive, because 
of the instabilities that occur as large time step sizes are used [8]. 
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When the swapping of edges occurs, the solution, which is stored at the vertices of the 
triangulation, needs to be modified. For capturing weak solutions of hyperbolic conserva- 
tion laws, conservation in time is a requirement. If conservation were not an issue, the 
solution need not be modified since it merely changes from one piecewise-linear representa- 
tion to another and the two are acceptable second-order accurate solutions. Note that the 
swapping takes place at a fixed time step after the solution has already been evolved up 
to that time level. It is possible to carry out another iterative loop so that the unsteady 
residual is driven to zero in the new configuration. This would double the work done at 
each timf step. Instead, we propose a noniterativq conservative interpolation at each time 
step. With the lumped mass matrix, Ecpi. (4), the requirement for conservation is that 
Vertices WV be conserved before and after the swapping. When an edge is swapped, the 
control volumes change for all the four points forming the quadrilateral. Figure 15 shows 
a four point quadrilateral subset of a triangulation. The initial triangulation given by the 
triangles 125 and 126 and the triangulation after swapping given by the triangles 156 and 
256. The portions of the control volumes associated with the four vertices that fall inside the 
quadrilateral region are illustrated as well. A conservative interpolation can be performed 
by treating the solution inside the control volumes as piecewise constants and computing 
geometrically the fractions of the old control volumes comprising the new ones. This would 
require complex intersection computations, especially for the vertex-centered scheme,, where 
the control volume edges are segmented edges. This algorithm may be viewed as a particular 
application of the algorithm due to Ramshaw [32] . An algebraic approach that only involves 
areas A\, A 2 , A 5 , A 6 , Ay, Ay, Ay, Ay is attractive, but the system becomes underdetermined 
if conservation is assumed. The drawbacks of these algorithms are algorithmic complexity 
and their diffusive character. The latter is particularly nettlesome because it degrades the 
accuracy to first order. For example, swapping back to the initial configuration and repeating 
the procedure in Figure 15, will result in equal values at the four vertices. 

5 5 


1 


6 6 

OLD NEW 

Figure 3: Two possible triangulations of four points and control volumes of vertices. 
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In the following, a conservative, linearity-preserving interpolation procedure is presented. 
The derivation is given assuming that the lumped mass approximation is adopted. We 
observe that the following identity holds: 

WV= Wv ’ ( 35 ) 

Vertices Triangles 

where V is the area of the control volume and V is the area of the triangle, and W, the 
average value in a triangle, is given by the average of the vertex values. Referring to Figure 
15, when an edge is swapped, the contributions from the triangles, 125 and 126, change to 
the contributions from the two swapped triangles, 156 and 256, on the right hand side of 
Eqn. (35). We further note that each term in the right hand-side may be viewed as volume 
under a “roof” in the x — y — W space. If the data were linear, the volumes under the two 
triangulations would be identical. In that case, no changes need to be made to the solution 
at the vertices. In the general case, we compute the contributions to the right-hand side of 
Eqn. (35) from the two triangulations as: 

T 0 i d = ^ [(Wi + W 2 + 1E 5 )V 125 + (lEi + W 2 + 1E 6 )V 126 ] , 

Tnew = g [(hEl + IE 5 + hEe)Vl56 + (IE 2 + W 5 + !T6)V256] • (36) 

Changes are now made to the vertex values in a conservative manner by distributing the 
difference, T 0 u — T neW} to the nodes. In principle, the changes can be made to all the four 
vertex values and would still result in a linearity-preserving, conservative scheme, but in 
inspecting Figure 15 we notice that with swapping, the control volumes associated with 
vertices 1 and 2 (connected by the old edge) can only decrease, whereas those associated 
with vertices 5 and 6 (connected by the new edge) can only increase. Therefore changes need 
to be made only to vertices 5 and 6. We apportion the difference, T 0 u — T neW} equally: 

jjrnewTznew jjr t inew , ^ / rji rji \ 

w 5 v 5 = W 5 v 5 + -{J-old — J-new), 

W ne W yne W = ^(T old - T new ). (37) 

We now show that the interpolation formulas given by Eqn. (37) satisfy the conservation 
property. First note that the overlapping control volume ry, given by the union of the 
triangles meeting at vertex z, is related to the nonoverlapping control volume R by the 
formula, 

v t = 3V t = J2V. (38) 

3 

where the sum is over all the triangles meeting at vertex i. After swapping, we have 

W 5 new = v n 5 ew = v 5 -V 125 + V 15 e + V 25 e = v 5 + V 126 . (39) 

We can thus derive the following relations: 

v 5 new = V + V 126 /3. 
v 6 new = V 6 + V 125 /3, 

vr w = U-V256/3, 

v 2 new = v 2 - V 156 /3. 
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With the help of these relations, it is easy to show that the conservation property holds: 


yynewynew + yynewynew + ^ynew 


+ w 2 v 2 new 


W 5 V 5 + W 6 V 6 + W 1 V 1 + W 2 V 2 , 


( 40 ) 


where use has been made of the fact that the values at vertices 1 and 2 do not change. 

In the formulas given by Eqn. (37), the difference T 0 u — T new has been apportioned equally 
to the two vertices. It is easy to see that this could introduce new extrema in the solution. 
The degree of freedom available in how this apportionment is made to the two vertices can 
be used to effectively prevent new extrema. Let 


W max = Max(W u W 2 ,W 5 ,We) 
W mm = Min^^IL^We). 


For i = 5 or 6 compute, 


n = 


V r 

r l 

V? 


? w max -w ■ 

Tnld~T new 

! 

T n id~T new 


if T 0 i d ~ T new > 0 
if T 0 i d ~ T new < 0. 


Let r = Min[0.5, r 5 , r 6 ]. The interpolation formulas can be written in a compact form as 
follows: 


W new 

W new 


W i T nld -T r , 
» » 5 — r ynew 
v 5 

T/fy _i_ Toid-T, 

6 — f" t yneu 


(r-r 5 )(r-r 6 ) r, r , (r-0.5)(r-r 6 ) , (r-0.5)(r-r 5 ) _ \ 

(0.5-r 5 )(0.5-r 6 ) U ' ' (i’5-0.5)(r 5 -r 6 ) ' (r 6 -0.5)(r 6 -r 5 ) V /_ 

(r-r 5 )(r-r 6 ) r, r , (r-0.5)(r-r 6 ) r, \ , (r-0.5)(r-r 5 ) ' 

(0.5-r 5 )(0.5-r 6 ) U ' ' (r 5 -0.5)(r 5 -r 6 ) V / ' (r 6 -0.5)(r 6 -r 5 ) 


. ( 41 ) 


The formulas are conservative, linearity-preserving, and also do not introduce new extrema. 
In contrast to the conservative formulas that assume piecewise constant values and result in 
equal values at the vertices upon repeated application, the new formulas converge locally to 
a linear (possibly least squares) profile for the four data points upon repeated application. 
The swapping and the interpolation formulas generalize to three-dimensional tetrahedral 
tessellations and are discussed in the Appendix. 


7 Results 


First, results from a one-dimensional example are presented illustrating the role of the mass 
matrix. On a uniform mesh, since the vertex and the centroid of its control volume coincide, 
the mass matrix can be lumped without suffering any adverse consequences. The situation is 
different if a mesh with variable mesh widths is considered. The one-dimensional advection 
equation 


du du 
dt dx ’ 


( 42 ) 


is solved on a random grid. The scheme stores the pointwise values U{ at locations aq. The 
initial condition is a Gaussian and the profile is advected by marching to a fixed time. A grid 
refinement study is carried out using a constant CFL number of 0.3. The spatial derivative 
is approximated in a MUSCL scheme [40] as 


_ u i + 1/2 U i- 1/2 

y dx J . Ax 
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where uf +1 j 2 , the value interpolated to the left side of the face i + 1/2, is given by 


U f+l/2 = Vi + (x i+1 /2 ~ Xi) U% + 1 _ U% 1 ■ (44) 

%i -\- 1 %i — 1 

The mass matrix, which is tridiagonal, is inverted using the Thomas algorithm. We have 
experimented with two definitions of the mass matrix. The first one derives the mass matrix 
by assuming a piecewise linear distribution of data between the grid points and computes 
the average over the control volume [x 8 _i/ 2 , ay +1 / 2 ]: 

1 ay - ay_i 3 1 x i+1 - ay 

T u i - 1 + 7 u i + 7 u i + 1 

4 x i+1 - ay_i 4 4 x i+1 - ay_i 

A second definition of the mass matrix is derived by computing the average of the recon- 
struction polynomial within a control volume. This polynomial is given by 



u(x) = Ui + (x — Xi) 

The average over the control volume is given by 


^2 + 1 1 
X%-\-\ Xi — i 


Xi - 1 — 2 ay + ay+i ay_i — 2 ay + ay+i 

-n 8 _l + Ui -| — ; r n 8 _i 


8(ay+i + ®i-i) 


8(®i+i + ®i-i) 


(46) 


(47) 


Figure 1 compares the errors in L 2 norm with the mass matrices given by Eqn. (45) and 
Eqn. (47), and with the lumped mass matrix. All the schemes exhibit second order accuracy 
but the errors are larger with the mass matrix given by Eqn. (45). The results obtained 
with the lumped mass matrix are almost identical to those obtained with Eqn. (47). Taylor 
series expansion would imply a first order error with the lumped mass matrix on a random 
grid, whereas Figure 1 clearly indicates second order accuracy. The results therefore reveal 
the inadequacy of local analysis. The results obtained with the usual finite element mass 
matrix, with hi = ay — ay_i, 


2 hj 

6 (hi + hi + 1) 


4 

u i - 1 + ~7. u i + 

6 


6 (hi + hi + 1) 




( 48 ) 


are also shown in Figure 1 and again display larger errors compared to the lumped mass 
approximation. The reason for this is that the finite element mass matrix is consistent 
with a Galerkin method which corresponds to a central difference discretization, whereas 
the spatial differencing employed here is upwind-biased. After experimenting with a one- 
parameter family of mass matrices, we have found that the lumped mass matrix gives the 
lowest errors with this particular spatial discretization. 

It is well known in finite element literature [38] that in some cases the lumping of the 
mass matrix does not compromise the solution accuracy, but that the mass matrix may play 
a crucial role when higher-order discretizations are considered. To examine this, we employ 
a quadratic reconstruction procedure utilizing point values. With hi = ay — ay_i, we obtain, 


U i+ 1/2 — 


+ 


h 


2 

4 + 1 


4 hi(hi + hi. |_i) 


Ui-l 


1 . h ^ |-i \ hi + 1 + 2 hi 

~ 4 — —j — u i + Try ; — r~: u i- 

2 4 hi ) 4(/q + i + hi) 
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The finite volume mass matrix is derived by determining the average of the quadratic dis- 
tribution and is given by: 

ctiUi-i + a2iii + a^Ui + i, 

_ — 2/i?_|_i + 2/q/q_|_i + lij 

12 hi(hi + /q+i) 

_ hf-|_ i + 2/q/q_|_i — 2 lij 

12/q+i (hi + /q+i) 

a 2 = 1 — oq — a 3 . 

The standard Runge-Kutta scheme which is fourth order accurate in time is used for the 
higher order computations. Figure 2 shows the error plots using the lumped mass matrix 
and the full mass matrix on the random grid. It shows that with the lumped mass matrix, 
the accuracy eventually degrades to second order as the grid is refined, whereas using the 
full matrix yields the third order accuracy of the spatial discretization. We have observed 
that using any other definition for the mass matrix degrades the accuracy to second order. 

The implications for the scheme in multiple dimensions are clear. As long as only a 
second order accurate scheme is used and we operate with either cell- vertex or cell-centered 
data, the mass matrix may be lumped without any loss of order of accuracy. The mass 
matrix can also be ignored for second (and higher) order accurate schemes if a strict cell- 
average interpretation is employed. If point values are used to construct third and higher 
order accurate schemes, the accuracy will degrade if the mass matrix is lumped. For higher 
order accurate schemes based on point values, the indirect mass matrix inversion technique 
discussed earlier will help preserve the order of accuracy of the scheme. 

We next present results from two-dimensional inviscid calculation over a pitching airfoil. 
The transonic flow is over a sinusoidally oscillating NACA0012 airfoil where the angle of 
attack ait) varies according to the formula 

ait) = a m + a 0 sin(iot) (49) 

For the test case chosen, a m = 0.016°, a 0 = 2.51°, k = = 0.0814 and the freestream 

Mach number, M = 0.755. Computing this flow using an explicit scheme is time-consuming 
because of the low frequency. The flows are computed using two meshes, referred to as GRID1 
and GRID2, each having 6336 vertices. These are shown in Figures 3 and 4, respectively. 
GRID1 is generated by drawing diagonals in a structured C-mesh and is fairly uniform. 
GRID2 is generated by random perturbations on GRID1. In GRID2, the centroids of the 
control volumes formed by the median dual are not represented well by the vertices. Figure 5 
shows the lift histories during the third cycle of oscillation. Four curves are shown, namely, 
the histories with the lumped and full mass matrices for GRID1 and GRID2. Since a 
projection-evolution approach is not followed, the mass matrix is derived by using a definition 
similar to Eqn. (45). As expected, the mass matrix has little impact on the integrated 
quantities even in the random mesh. The differences in the solutions between the two grids 
are likewise insignificant. The CPU time increases by about 15% when the full mass matrix 
is included. These examples have been run with a maximum physical CFL number of 500, 
corresponding to using 54 time steps per sinusoidal oscillation of the airfoil. The number of 
iterations for the inner multigrid procedure is fixed at 30. Figure 6 shows the convergence 
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Figure 6: GRID1 about an NACA0012 airfoil with 8338 vertices. 



Figure 7: GRID2 about an NACA0012 airfoil with 6336 vertices. 
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of 4-grid agglomeration multigrid procedure during a particular time step with thq lumped 
and the full matrices where the L- 2 norm of the unsteady residual R* is plotted as a function 
of the multigrid cycles using 4 grid levels. The convergence improves slightly when the mass 
matrix is included. The reason for this is that the mass matrix can be recast as an implicit 
smoothing operator. The remaining examples have been computed with the lumped mass 
matrix. Finally, Figures 7 and 8 shows the effect of the physical time step size. Three 



Alpha 

Figure 8: Lift histories during the third cycle of motion. 

lift and moment histories are shown, employing 40, 20 and 10 steps per pitch cycle using 
the lumped mass matrix and deforming grids. The grid velocities arq computed by using 
Ecpi. (27). The nodes are repositioned by the procedure described earlier. No swapping 
of edges occurred in this computation. The number of multigrid cycles used at each time 
step are respectively, 15, 20 and 20. With 10 steps per pitch cycle, somf ■discrepancy may 
be observed; increasing the number of multigrid cycles (thus solving the nonlinear problem 
better at each time step) did improve the solution. Also shown is a comparison with the 
experimental data of Landon [16] as well as with a structured grid computation employing 
a TVD scheme on the same grid [42], The differences in the moment histories between the 
structured and unstructured grid computations may be attributed to the differences in the 
spatial discretization. To assess the effects of swapping, we modify the Delaunay criterion to 
force the swapping of edges. Given a quadrilateral with two possible triangulations, rather 
than swap to the new configuration if it has a larger minimum angle, we swap if some 
multiple (1.1 in this example) of the minimum angle in the new configuration is larger than 
the minimum angle in the old configuration; only one pass of this algorithm is performed. 
On a coarse triangular grid generated from a structured 128 X 32 grid, the moment histories 
during the first three pitch cycles without the swapping are shown in Figure 9. Figure 
10 shows the histories during the first five pitch cycles with the swapping of edges. It is 
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seen that the transient response is indeed quite different especially during the second cycle, 
although the response does become periodic during the fifth cycle. Figure 11 shows the 
moment histories on a triangular mesh generated from a 256 X 64 grid with the swapping of 
edges. The transient response is not nearly as erratic as on the 128 X 32 grid with swapping 
and more in line with the response on the coarse grid without swapping. Note that the 
swapping and subsequent conservative interpolation introduces errors proportional to Ax 2 , 
which decrease as the grid is refined. 

We next present a laminar unsteady calculation of impulsive start-up flow over a cylinder 
at a Reynolds number of 1200. Rumsey [34] has computed this flow using a structured grid 
code and has made detailed comparisons with experimental data [26]. We employ the same 
grid (192 X 64) as in [34], but divide the quadrilaterals into triangles. The flow separates 
and eventually vortices are shed. The computed Strouhal number is 0.225 as compared to 
the experimental value of 0.215 and the computed value of 0.222 using the structured grid. 
The time is nondimensionalized as i = t/(d/Uoo) where d is the diameter of the cylinder 
and Uoo is the freestream velocity. The sequence of non-dimensional time steps At is chosen 
as in [34]: At = 0.01 for t up to 6.0, At = 0.02 between t = 6.0 and 9.0 and At = 0.05 
above i = 9.0. We use the agglomeration multigrid strategy with six grids where the edge- 
coefficients needed for computing the viscous and inviscid terms on the coarse grids are 
precomputed as in [22], Seven multigrid iterations were used at each time step yielding 
about of 2-3 orders of reduction in the unsteady residual. The centerline velocity at a time 
i = 2.9 is plotted in Figure 8, and the maximum reverse flow velocity is plotted in Figure 9 
as a function of i during the initial bubble growth phase. Good agreement may be observed 
with experiment. There is some discrepancy with the structured grid code which may be 
on account of the fact that it did not employ inner iterations. Thus, errors arising from 
factorization, linearization and and mismatch of operators may be present in the structured 
grid solution. 

The final test case represents an exploratory inviscid computation of flow over a multi- 
element airfoil system undergoing deployment, as an example of bodies in relative motion. 
The geometry represents a sectional cut of the wing of the NASA Langley Transport System 
Research Vehicle (TSRV), and was obtained by direct measurement of the full scale aircraft 
(Boeing 737-100) [46]. The initial position, depicted in Figure 10, corresponds to the 15° 
flap setting, while the final condition represents the 40° flap setting, which is the approach 
configuration for this high-lift system. The initial grid about the 5-element airfoil is generated 
using the advancing front Delaunay triangulation method [21] at the 15° flap setting and is 
displayed in Figure 10. The grid has 26,191 vertices. Figure 11 displays the Mach contours 
for steady flow at this setting at an angle of attack of 5° and a freestream Mach number of 
0 . 2 . 

We compute the time-accurate flow solution to full deployment using a five-grid agglomer- 
ation multigrid procedure. The motion is prescribed as a linear variation of both translation 
and rotation of the various airfoil elements, and is computed in 300 time steps. The non- 
dimensional time for deployment, defined as t = // \Jpoo / Poo , is taken as 10. Here / is the 
chord of the main element and p ^ , p ^ are the freestream pressure and density. The grid 
restructuring involved presmoothing, spring analogy and edge swapping. Figures 12 and 13 
depict a closeup of the grid in the flap region and the instantaneous flow solution at full de- 
ployment. Grid quality, while not as good as with the initial mesh, is acceptable considering 
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the large scale motion. For this case, rather than the 3-point backward approximation, we 
have chosen to use the implicit midpoint (trapezoidal) rule: 


yn+l W n+l _ yn W n R(]y n+1 ) + R(W n ) 

At + 2 “ 

When the connectivity of the mesh changes, the solution changes from one piecewise linear 
representation to another and the new control volumes are used to evolve the solution to 
the next time step. Using a 3-point backward rule would require that the volumes at the 
previous time step be changed as well. Figure 14 plots the total lift history as well as the lift 
histories of the elements as a function of physical time. The total lift, based on the chord 
length in the fully nested position, increases from an initial value of 2.53 to a final value of 
3.10. 


8 Conclusions 

An efficient implicit time integration procedure has been developed. The implicit system is 
solved by using the agglomeration multigrid procedure. The issue of the mass matrix which 
arises in cell- vertex methods is addressed. It is shown that lumping of the mass matrix may 
be done for second-order accurate schemes without any degradation in order of accuracy. For 
higher order methods based on point values, it is shown that the lumping of the mass matrix 
degrades the order of accuracy of the scheme. Inviscid and viscous calculations have been 
presented for flow over a pitching airfoil and an impulsively started cylinder and the results 
have been compared with experiments and other computations. A mesh point movement 
strategy has been proposed and tested. This has been used to compute inviscid flow over a 
multi-element airfoil system undergoing deployment. 
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Figure 11: Moment histories during the third cycle of motion with with 40, 20 and 10 steps 
per cycle. 
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Figure 12: Moment histories on the coarse O-mesh without swapping of edges. 
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Figure 13: Moment histories on the coarse O-mesh with forced swapping of edg 
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Figure 14: Moment histories on the fine O-mesh with forced swapping of edge; 





0.5 


Unstructured grid 

Structured grid [29] 



Figure 15: Velocity distribution at t = 2.9 on the symmetric axis of the wake. 



Figure 16: Time history of the maximum reverse flow velocity on the symmetric axis of the 
wake. 
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Figure 17: Initial grid at the semi-deployed position (15° flap deflection). 



Figure 18: Mach contours for steady flow at the initial position (15° flap deflection). 
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Figure 19: Closeup view of the grid at full deployment (40° flap deflection). 



Figure 20: Mach contours of the instantaneous solution at full deployment (40° flap deflec- 
tion ) . 
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Appendix 

We discuss the swapping and interpolation procedures in three dimensions. Given 5 points 
in three space, the region enclosed by the convex hull can be tetrahedralized in one of three 
ways: 

1. 4 tetrahedra sharing a common interior vertex. 

2. 2 tetrahedra sharing a common triangular face. 

3. 3 tetrahedra sharing a common edge. 

Swapping in three dimensions consists of switching between the second and the third choices 
to improve the quality of the tetrahedralization [17]. Figure 25 illustrates the 2- and 3- 
tetrahedra configurations . 
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(a) 2-tetrahedra configuration (b) 3-tetrahedra configuration 

Figure 22: Two possible tetrahedralizations of five vertices. 

Conservative, linearity-preserving interpolation, when switching from the 2-tetrahedra 
configuration to the 3-tetrahedra configuration, is similar to the two-dimensional situation, 
since only the values at the vertices joined by the new edge need to be changed. Define 

T 2 = i [( W 1 + W 2 + W 3 + W 4 ) V 1234 + ( W 1 + W 2 + W 3 + W 5 ) v 1245 ] , 

T 3 = — [( W\ + W -2 + W 4 + MG ) V1245 + (IT2 + W 3 + W 4 + MG ) V2345 + (IT3 + W\ + W 4 + MG ) V314.5] . 

Changes are now made to the vertices 4 and 5 in a conservative manner by distributing the 
difference T - 2 — T 3 using formulas similar to Ecpi. (41), where the W max and W m - m are taken 
as the maximum and minimum over the five points. 
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Swapping from the 3- to the 2-tetrahedra configuration is different, however. Here changes 
are made to the vertices 1,2 and 3 by distributing the difference T 3 — T 2 . Thus, there are 
two degrees of freedom, which are needed to enforce the condition that no new extrema be 
created. For z = 1,2,3 compute 


Vi = 


> W max -w, 


T /new v v max 

Ts-T : 


vr 


•W„ 


-h 


t 3 -t 2 ’ 


if T 3 
if T 3 


T 2 > 0 
T 2 < 0. 


and define 

r = Min[ri,r 2 ,r 3 , -]. 

Now assume r = r\. Then 


Dehne s - 

W new 

w 3 ew 


jTrnewj rnew jit -\rnew , / rp rji \ /-rrnew 

W 1 V 1 = W 1 V 1 + r(l 3 — 1 2 )/V 1 
Min[r 2 , r 3 , 0.5] . Changes are made to vertices 2 and 3 as follows: 

■ Wo 4- Ts ~ T2 n — r) l ( s ~ r H(s-’-3) n c , (s — 1 0.5)(s — 7-3 ) , (s — 0.5)(s — r 2 ) _ \1 

2 W ynew t ) |_ (0.5 -r 2 ) (0.5 -r 3 ) U ■ ' (i’2-0.5)(r 2 -r 3 ) ' (r 3 -0.5)(r 3 -r 2 ) V )\ ’ 

- Wo 4- T 3~ T2 (1 — r) l ( s ~ r2 )( s ~ r 3) n c I (g — 0.5)(s — r 3 ) ^ \ I (s-0.5)(s-r 2 ) 

3 ' ynew t ) |_ (o.5-r 2 )(0.5-r 3 ) U ' ' (r 2 -0.5)(r 2 -r 3 ) V / ' (r 3 -0.5)(r 3 -r 2 ) 


Formulas for the remaining cases, when r = |,r 2 ,r 3 , can be derived in a similar man- 
ner. 
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